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We have developed a numerical differentiation scheme which eliminates evaluation of overlap determinants in 
calculating the time-derivative nonadiabatic couplings (TDNACs). Evaluation of these determinants was the 
bottleneck in previous implementations of mixed quantum-classical methods using numerical differentiation 
of electronic wave functions in the Slater-determinant representation. The central idea of our approach is, 
first, to reduce the analytic time derivatives of Slater determinants to time derivatives of molecular orbitals, 
and then to apply a finite-difference formula. Benchmark calculations prove the efficiency of the proposed 
scheme showing impressive several-order-of-magnitude speedups of the TDNAC calculation step for midsize 
molecules. 


Simple quantum-classical methods for simulating nona¬ 
diabatic dynamics, Ehrenfest and fewest-switches surface 
hopping (FSSH), 1,2 often provide accurate and efficient 
ways to investigate chemical processes involving several 
electronic states. The simplicity of these methods stems 
from restricting quantum mechanical consideration to the 
electronic part and treating the nuclear part classically 
with minimal intervention of quantum mechanics. For 
describing the quantum evolution of the electronic sub¬ 
system its non-stationary electronic wave function is writ¬ 
ten in terms of the adiabatic eigenfunctions {Tj(r;R)} 
of the electronic Hamiltonian H e as 

i>{t, r; R(i)) = X c -/(*)'M r ; R W)- (!) 

j 

The time-dependent coefficients cj{t) then can be ob¬ 
tained via projecting the time-dependent electronic 
Schrodinger equation onto the orthonormal basis of 
{Tj(r;R)} (atomic units are assumed hereinafter) 

3-^- = X cj(5kjEj(R) - irpcj ), (2) 

where 5 kj is the Kronecker delta, Ej( R) are the adia¬ 
batic potential energy surfaces (PESs), and 

TKJ = {^K\dt*j), K^J, (3) 

are the time-derivative nonadiabatic couplings (TD¬ 
NACs). Using the chain rule one can further decompose 
t kj = R • d K j, where d KJ = (>P K | VrTj) is the 3M- 
dimensional (M is the number of nuclei in the system) 
derivative couplings vector, and R is the 3M-dimensional 
nuclear velocity vector, d kj are implemented analyti¬ 
cally for some electronic structure methods, such as multi- 
configurational self-consistent field (MCSCF), 3,4 config¬ 
uration interaction singles (CIS), 5,6 and multi-reference 
configuration interaction (MR-CI).'’ 8 However, many 
electronic structure methods either lack of the analytic 


implementation ( e.g. XMCQDPT 9 ) or have intrinsic 
problems in their definition. 10-15 Of course, numerical dif¬ 
ferentiation is always an option for evaluation of d^j’s 
but it is also quite computationally expensive considering 
the dimensionality of these quantities. 

On the other hand, it has been recognized that it is 
more efficient to apply numerical differentiation to TD¬ 
NACs directly. 16 For example, any of the following for¬ 
mulae can be employed 

Tkj =-^{'S’K( t )\'S>j(t +At)) + o(At), (4) 

TKJ = JKt ( ^ I + At ^ 

- {V K (t + At) | j(t - At))) + o(At 2 ), (5) 

giving rise to the first-order forward and the second-order 
central finite difference schemes, respectively. Moreover, 
as was shown in Ref. 17, use of numerical TDNACs can 
be advantageous close to surface crossings. However, in 
the conventional FSSH method, d kj quantities are also 
needed to rescale nuclear velocities if a hop between elec¬ 
tronic surfaces Ek{ R) and -E’j(R) occurs. 

Recently, to avoid numerical evaluation of Akj , a sim¬ 
pler version of the FSSH method has been suggested. 18 In 
this simplified version nuclear velocities are rescaled uni¬ 
formly after a surface hop. It was shown that this simpli¬ 
fied scheme can adequately model nonadiabatic dynamics 
and deviates from the regular FSSH algorithm only in re¬ 
gions where dA'j’s change rapidly, but these deviations 
have only a minor effect on overall dynamics. 19 Thus, if 
we focus on the simplified FSSH method, the only re¬ 
quired nonadiabatic coupling terms will be TDNACs. 

In commonly used numerical formulations [Eq. (4) or 
(5)], TDNACs are obtained from an electronic overlap 
matrix 

£ Kj{t , t + At) = (T k( i) | ’P j(t + At )), (6) 

whose matrix elements require Slater-determinant pair 
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overlaps for electronic wave functions in the Slater- 
determinant representation 19 

^Kj(t,t + At) = E 

{p,q} 

x {®{ P }(t)\®{ q }(t + At)}, (7) 

where Cy p y and Cfa are coefficients of the |3>{ p }(i)) 
and |${ g }(t + At)} Slater determinants. Here, we use 
collective indices {p} and {< 7 } denoting sets of or¬ 
bitals present in the Slater determinants | < E > { p }(t)) and 
|d>{ 9 }(t + At)). This scheme quickly becomes computa¬ 
tionally expensive with the system size, because it re¬ 
quires evaluation of many Slater-determinant pair over¬ 
laps ( < f»{ P }(t)I*&{(}} (t + At)) given by the Lowdin for¬ 
mula 20,21 

($ {p }(t)|$ w (t + At)) = det S[{p}{ 9 }], ( 8 ) 

where S[{p}{< 7 }] is the overlap matrix of molecu¬ 
lar orbitals comprising the determinants |d > { p }(t)) and 
|$ M (t + At)). The Lowdin formula appears as a result 
of non-orthogonality between sets of orbitals at different 
times. The computational cost of detS[{p}{q}] calcula¬ 
tion grows cubically with the number of electrons in the 
system, N e . 22 Considering that all pairs of determinants 
in Eq. (7) need to be evaluated, use of Eq. (7) in large 
systems makes the evaluation of TDNACs the bottleneck 
of mixed quantum-classical simulations. 

In this Letter we show how computing of the determi¬ 
nant overlaps can be avoided in numerical evaluation of 
TDNACs without introducing any approximations and 
by making the procedure faster by at least a factor of 
~ A r ° cc for each determinant pair (N occ = N e /2 for the 
closed-shell case). We illustrate the performance of our 
approach by computing TDNACs at the CIS level of the¬ 
ory, which is one of the simplest methods for treating ex¬ 
cited states. Our developments can be straightforwardly 
applied to any other method presenting wave functions as 
linear combinations of Slater determinants (e.g., MR-CI 
or MCSCF). 

In the CIS method, excited-state wave functions are 
written as linear combinations of singly-excited Slater 
determinants |<I>“) = d^di |$ 0 ) obtained from the ground- 
state Hartree-Fock determinant |$o) 

l^) = E C ^l $ “)> (9) 

ia 

with coefficients C(£ defined by the secular matrix prob¬ 
lem H e C = EC. Here, we follow the common conven¬ 
tion where subscripts a, 6, c,... denote virtual orbitals, 
i,j, k,... label occupied orbitals, and p,q,r ... are used 
for any type of orbitals. 

To avoid computational difficulties associated with 
overlap determinants [Eq. ( 8 )] we will start with the for¬ 
mal definition of TDNACs as time derivatives [Eq. (3)] 
and postpone applying a finite difference scheme until we 


account for the anti-symmetric structure of Slater deter¬ 
minants. Assuming real-valued molecular orbitals and 
CIS coefficients, TDNACs can be written as 

TKJ = E(OtC/ 6 <$“ | $*> + c£ Cj b <$“ | d t fa>)) ■ 

iiab 

( 10 ) 

All terms in Eq. (10) refer to the same t, hence, there 
are no complications with orbital non-orthogonality as 
in Eq. (7). One may apply the Slater-Condon rules to 
the first term, but not to the second one, since the time 
derivative dt is not an operator in the space of electronic 
variables. Instead, we differentiate determinants |$*) di¬ 
rectly 

d t |d>5> = E !<> + !<>, (H) 

where the notation ) means that a molecular orbital 
t f> p is replaced with the time derivative dfaq. Therefore, 
TDNACs between determinants become 

I dt $$> = E <*“l*Sfc> + (*“!<> • ( 12 ) 

The last term in Eq. (12) is reduced to Sij fa a \dtfa), 
while only one term with k = i and a = b from the 
sum over k survives because of orthogonality conditions: 
{(j) p \dt(j)p) = 0 (for real orbitals) and (fa\4> q ) = 5 pq . Fi¬ 
nally, we have: 

($? | &$$) = S i:i {ct> a | dtfa) - P l3 S ab (fa | d t fa) , (13) 

where Pij is an additional phase factor which depends on 
the ordering convention for the orbitals in the Slater de¬ 
terminants. There are two common choices which result 
in different Pij values: 

\$i) =det{...,fa- 1 ,fa,fa +1 ,...}, P^ = 1, (14a) 

l$“) = det {..., fa-ufa+i ,..., fa}, Pij = (-l) 1 ^, 

(14b) 

Use of Eq. (13) leads to substantial reduction of 
TDNAC computation scaling because Eq. (10) is simpli¬ 
fied to 

TKJ = E cf a d t c( a + E C?C& (fa I dtfa) 

ia iab 

-Y,PnC£c] a (fa \d t fa). (15) 

ija 

Each term in Eq. (15) scales as N occ N v i rt , N occ N 2 irt , 
and A 2 cc lVvirt, respectively, where N occ and N v j rt are the 
numbers of occupied and virtual orbitals. In typical cal¬ 
culations, N v j rt > N occ , and the second term [Eq. (15)] 
is dominating in the overall computational cost provid¬ 
ing the overall cubic scaling with the size of the system. 
This scaling should be compared with the N^ cc N 2 1It scal¬ 
ing of Eq. (7) for the CIS method with the determinant 
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scheme. Note that evaluation of all terms in Eq. (15) can 
be reformulated as highly efficient matrix-matrix multi¬ 
plications. 

The possibility of reducing time derivatives of deter¬ 
minants to time derivatives of orbitals has been already 
mentioned in several works. 23,24 However, Eq. (13) has 
never been derived explicitly; our treatment, therefore, 
provides a rigorous foundation for the orbital formula¬ 
tion and for its extension to multi Slater determinant 
wave functions. 

To apply Eq. (13), one has to convert it into a 
corresponding finite-difference expression. Any finite- 
difference expression requires continuity of orbitals at dif¬ 
ferent times. However, orbital phases at different times 
are arbitrary , reflecting the existence of the wave func¬ 
tion gauge (phase) degree of freedom. Thus an appro¬ 
priate orbital phase matching and tracking procedure is 
necessary. 

The finite-difference counterpart of Eq. (13) is obtained 
by substituting 

(</)p\dt(j) q ) -A -^S pq {t,t + At), p^q (16) 

where 

S pq (t,t + At) = {4> p {t)\4> q {t + At)) (17) 

is the orbital overlap matrix. To keep track of relative 
signs of orbitals at t and t + At we introduce an integer¬ 
valued matrix O, which is obtained from S (£, f + At) by 
rounding off its matrix elements to ±1 or 0. O has a 
structure of the signed permutation matrix as long as At 
is sufficiently small. 20 Performing the permutation and 
sign changes of molecular orbital coefficients in C(f +At) 
according to O we obtain matrix C(f + At). The set 
of orbitals C(t + At) is subsequently used to calculate 
the CIS coefficients at the moment t + At. Tracking and 
phase matching for the CIS states remain the same as for 
the determinant-based procedure. 19 Computational over¬ 
head for the orbital tracking and phase matching scales 
as (IVocc + Avirt) 2 and is negligible compare to other com¬ 
ponents of the algorithm. 

To test the accuracy and efficiency of the proposed 
scheme, we benchmark it against the conventional scheme 
based on Eqs. (7) and (8) as implemented by Pittner and 
coworkers in the Newton-X program. 26,27 Table I shows 
that accuracies of the orbital- and determinant-based 
schemes are very similar as could be expected from nu¬ 
merical schemes of the same order. For the efficiency com¬ 
parison it is worth noting that the Newton-X scheme 
has been used with a screening threshold of 5 x 10 -5 for 
the products of CIS coefficients to reduce the number 
of determinant overlaps in Eq. (7). Our implementation 
uses matrix-matrix multiplication and thus do not em¬ 
ploy a screening procedure. Table II illustrates speedups 
achieved by the current scheme for two midsize organic 
molecules and various basis sets. The speedups are espe¬ 
cially prominent for small basis sets, where N occ is com¬ 
parable to the total number of basis functions. Increase 


TABLE I. Errors of numerical differentiation for representa¬ 
tive TDNACs in CH 2 NHj using the 3-21G basis. 28 TDNAC 
values converged up to 1 x 10 —7 , T31 = —1.170 x 10~ 3 and 
T21 = —2.375 x 10~ 4 are obtained with 0.05 fs time-step. 



Step size, fs 

Determinant 

Orbital 


0.1 

-1.4 x 10 -6 

-1.6 x 10 -6 

A731 

0.2 

-7.5 x 10" 6 

-8.1 x 10 -6 


0.5 

-5.3 x 10" 5 

-5.7 x 10 -5 


0.1 

-2.0 x 10 -7 

< 1.0 x 10 -7 

At 2 1 

0.2 

-7.0 x 10 -7 

-2.0 x 10 -7 


0.5 

-5.0 x 10 -6 

-2.5 x 10 -6 


of the basis set size makes the CIS coefficient product 
screening more productive, but the orbital-based scheme 
still outperforms the conventional scheme by more than 
two orders of magnitude. 

In practical applications one is more concerned with 
the total simulation time. Apart from TDNACs cal¬ 
culations, mixed quantum-classical nonadiabatic simula¬ 
tions include also the electronic-structure and classical- 
dynamics steps. For the electronic-structure CIS calcu¬ 
lations we have used the Gaussian program. 4 To give an 
idea of the overall speedup, we consider the first 50 fs of 
a single FSSH trajectory for the C 18 H 14 0 molecule us¬ 
ing the 6-31G** basis set and the 0.2 fs time-step. The 
electronic-structure part involved evaluation of character¬ 
istics for the three lowest electronic states. On a single 
Intel Xeon X5650 @ 2.67GHz CPU, it took 99 h in total 
to complete this trajectory for the original Newton-X 
procedure with 54 h 29 spent on the TDNAC calculations, 
the corresponding numbers for our algorithm are only 
45 h and 12 min. 

In conclusion, we have developed a numerical proce¬ 
dure which eliminates evaluation of overlap determinants 
in TDNACs calculations. . This elimination produces 
tremendous speedup in quantum-classical nonadiabatic 
simulations where evaluation of TDNACs was the bottle¬ 
neck. The central idea of our approach is to postpone 
introducing a finite-difference scheme, [Eqs. (4) or (5)] 
and to convert the expression for TDNACs given in terms 
of many-electron wave functions, Eq. (10), into the cor¬ 
responding orbital-based formula, Eq. (15). This alter¬ 
nation of steps allows us to manipulate with orthogonal 

TABLE II. Relative speedups for a single time step evaluation 
of TDNACs as compared to the Newton-X procedure. In all 
calculations TDNACs between the first three lowest electronic 
states were evaluated. 


Molecule 

Basis set 

Nocc 3, 

Xvirt 

Speedup 

C 18 H 14 O b 

STO-3G 

46 

44 

400 

6-31G** 

46 

290 

248 


cc-pVTZ 

46 

701 

172 

g 25 h 18 c 

STO-3G 

59 

59 

1372 


6-31G** 

59 

381 

292 


a Excluding Is core orbitals of C and O. 
b 4-(2-naphthylmethyl)-benzaldehyde. 
c 9- ((1-naphthyl)-methyl)-anthracene. 
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molecular orbitals and thus to avoid overlap determinants 
that arose as a result of orbital non-orthogonality in the 
determinant formulation. Benchmark calculations have 
proven the efficiency of the proposed scheme and illus¬ 
trated its potential for mixed quantum-classical studies 
of medium and large molecules. 
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